#include <iostream>
#include <vector>
#include <math.h>
//#include <lapacke.h>
#include "lapackSolver.h"

using namespace std;

struct Point // 定义点的结构存储坐标信息
{
    double x;
    double y;
};

class Grid // 网格
{
private:
    int N; // 划分
    int n;
public:
    vector<Point> G;
    Grid(int _N)
	{
	    N = _N;
	    n = (N+1)*(N+1);
	};
    Point get(int i, int j);
    vector<Point> get_allG(); // 所有网格点的坐标信息
    vector<Point> get_RDG(); // 求解器所需网格点
    vector<Point> get_RNG();
    vector<Point> get_IDG();
};

double Function0(Point g) // u的元函数
{
    double f;
    f = exp(g.y + sin(g.x));
    return f;
};

double Function1(Point g) // u的一阶导函数
{
    double f;
    f = exp(g.y + sin(g.x)) + cos(g.x)*exp(g.y + sin(g.x));
    return f;
};

double Function(Point g) // u的二阶导函数
{
    double f;
    f = exp(g.y + sin(g.x)) * (sin(g.x) * sin(g.x) + sin(g.x) - 2);
    return f;
};

class Value // 函数值
{
private:
    vector<Point> G;
    int N;
    int n;
public:
    Value(vector<Point> _G, int _N)
	{
	    G = _G;
	    N = _N;
	    n = (N+1)*(N+1);
      	};
    vector<double> get_allF0();
    vector<double> get_F0();
    vector<double> get_F1();
    vector<double> get_F2();
};

class Solver // 求解器
{
private:
    int N;
    vector<Point> G;
    vector<vector<double>> A_2;
    vector<double> F;
public:
    Solver(int _N, vector<Point> _G)
	{
	    N = _N;
	    G = _G;
	};
    vector<vector<double>> get_A();
    vector<double> get_F();
    vector<double> get_U(); // 解
    vector<double> get_U0(); // 解
    vector<double> yz_U(); // 验证
};
// --------------------Grid-----------------------------------------------------
Point Grid :: get(int i, int j)
{
    Point g;
    g.x = (double)i/N;
    g.y = (double)j/N;
    return g;
};
vector<Point> Grid :: get_allG()
{
    //vector<Point> G;
    G.resize((N+1)*(N+1));
    for (int i = 0; i < (N+1)*(N+1); i++)
    {
	G[i] = get(i/(N+1),i%(N+1));
    }
    return G;
};

vector<Point> Grid :: get_RDG()
{
    G.resize(n);
    for (int i = 0; i < n; i++)
    {
	G[i] = get(i/(N+1),i%(N+1));
    }
    return G;
};

vector<Point> Grid :: get_RNG()
{
    //vector<Point> G;
    G.resize((N+4)*(N));
    for (int i = 0; i < (N+4)*(N); i++)
    {
	if (i / (N-1) == 0)
	    G[i] = get(i+1,0);
	else if (i / (N - 1) == N+2)
	    G[i] = get((i-1) % N,N);
	else
	    G[i] = get((i-N+1)%(N+1),(i-N+1)/(N+1)+1);
    }
    return G;
};
//---------------------Value----------------------------------------------------
vector<double> Value :: get_F0()
{
    vector<double> F0;
    F0.resize(G.size());
    for (int i = 0; i < G.size(); i++)
    {
	F0[i] = Function0(G[i]);
    }
    return F0;
}
vector<double> Value :: get_F2()
{
    vector<double> F2;
    F2.resize(G.size());
    for (int i = 0; i < G.size(); i++)
    {
	F2[i] = Function(G[i]);
    }
    return F2;
}
// ---------------------Solver--------------------------------------------------
vector<vector<double>> Solver :: get_A()
{
    //vector<vector<double>> A_2;
    A_2.resize(G.size());
    for (int i = 0; i < G.size(); i++)
    {
	A_2[i].resize(G.size());
	if ((i / (N + 1) == 0) || (i / (N + 1) == N))
	    A_2[i][i] = 1;
	else if ((i % (N+1) == 0) || (i % (N+1) == N))
	    A_2[i][i] = 1;
	else
	{
	    for (int j = 0; j< G.size(); j++)
	    {
		if (i == j)
		    A_2[i][j] = 4*N*N;
		else if ((i == j+1) || (i+1 == j))
		    A_2[i][j] = -1*N*N;
		else if (i == j+N+1 || i+N+1 == j)
		    A_2[i][j] = -1*N*N;
		else
		    A_2[i][j] = 0;
	    }
	}
    }
    return A_2;
};

vector<double> Solver :: get_F()
{
    F.resize(G.size());
    for (int i = 0; i < G.size(); i++)
    {
	if (A_2[i][i] == 1)
	    F[i] = Function0(G[i]);
	else
	    F[i] = -Function(G[i]);
    }
    return F;
};

vector<double> Solver :: get_U()
{
    get_A();
    get_F();
    vector<double> U;
    U.resize(G.size());
    double *A1 = new double[G.size()*G.size()];
    double *F1 = new double[G.size()];
    for(int i = 0; i < G.size(); i++)
    {
	F1[i] = F[i];
	for (int j = 0; j < G.size(); j++)
	{
	    A1[i+j*(N+3)*(N-1)] = A_2[i][j];
	    cout << A_2[i][j] << "  ";
	}
	cout << "    " << F1[i];
	cout << endl;
    }
    
    lapack_int info,m,n,lda,ldb,nrhs;
    
    m = G.size();
    n = G.size();
    nrhs = 1;
    lda = G.size();
    ldb = G.size();
    
    info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m,n,nrhs,A1,lda,F1,ldb);
    
    for(int i = 0; i < U.size(); i++)
    {
	U[i] = F1[i];
    }
    return U;
};

vector<double> Solver :: yz_U()
{
    vector<double> U_0;
    U_0.resize((N+3)*(N-1));
    double *A1 = new double[(N+3)*(N-1)*(N+3)*(N-1)];
    double *F1 = new double[(N+3)*(N-1)];
    for(int i = 0; i < (N+3)*(N-1); i++)
    {
	F1[i] = F[i];
	for (int j = 0; j < (N+3)*(N-1); j++)
	{
	    A1[i+j*(N+3)*(N-1)] = A_2[i][j];
	    cout << A_2[i][j] << "  ";
	}
	cout << "    " << F1[i];
	cout << endl;
    }
    for (int i = 0; i < (N+3)*(N-1); i++)
    {
	U_0[i] = F1[i];
    };
    return U_0;
};

vector<double> Solver :: get_U0()
{
    vector<double> U0;
    get_A();
    get_F();
    U0 = F;
    double *A1 = new double[G.size()*G.size()];
    
    for(int i = 0; i < G.size(); i++)
    {
	for (int j = 0; j < G.size(); j++)
	{
	    A1[i+j*G.size()] = A_2[j][i];
	}
    }
    
    lapackSolver::solve(G.size(), G.size(), A1,& *U0.begin());
    return U0;
}; // 解
// -----------------------------------------------------------------------------
